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Abstract. In this paper we develop two conforming finite element methods for a fourth order 
bi-wavc equation arising as a simplified Ginzburg-Landau-type model for d-wave superconductors in 
absence of applied magnetic field. Unlike the biharmonic operator A^, the bi-wave operator is 
not an elliptic operator, so the energy space for the bi-wave equation is much larger than the energy 
space for the biharmonic equation. This then makes it possible to construct low order conforming 
finite elements for the bi-wave equation. However, the existence and construction of such finite 
elements strongly depends on the mesh. In the paper, we first characterize mesh conditions which 
allow and not allow construction of low order conforming finite elements for approximating the bi- 
wave equation. We then construct a cubic and a quartic conforming finite element. It is proved that 
both elements have the desired approximation properties, and give optimal order error estimates in 
the energy norm, suboptimal (and optimal in some cases) order error estimates in the and 
norm. Finally, numerical experiments are presented to guage the efficiency of the proposed finite 
element methods and to validate the theoretical error bounds. 
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1. Introduction. This paper concerns finite approximations of the following 
boundary value problem: 

(1.1) Sa^u-Au = f in f), 

(1.2) u = dnU ^ on dil, 

where < (5 ^ 1 is a given (small) number, 

Ou :— dxxU — dyyU, O'^u := □(□u), 

:= (rii, — rt2), dnU -.— Vu-n, 

fl C R^ is a bounded domain with piecewise smooth boundary dfl, and n :— (ni, 71,2) 
denotes the unit outward normal to dCl. As □ is the well-known (2-D) wave operator, 
we shall call the bi-wave operator throughout this paper. It is easy to verify that 

|--|2 / \ 9^"^ 2 

dx^ dxdy dy^ ' 

Hence, equation ( |1.1| is a fourth order PDE, which can be viewed as a singular 
perturbation of the Poisson equation by the bi-wave operator. As a comparison, we 
recall that the biharmonic operator is defined as 

AM-,y) :=A(Au(x,,)) = ^ + 2— -K^. 

Although there is only a sign difference in the mixed derivative term, the difference 

between and is fundamental because A^ is an elliptic operator while is a 
hyperbolic operator. 
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Superconductors are materials that have no resistance to the flow of electricity 
when the surrounding temperature is below some critical temperature. At the su- 
perconducting state, the electrons are believed to "team up pairwise" despite the 
fact that the electrons have negative charges and normally repel each other. The 
Ginzburg-Landau theory [9] has been well accepted as a good mean field theory for 
low (critical temperature) superconductors [11]. However, a theory to explain high 
Tc superconductivity still eludes modern physics. In spite of the lack of satisfactory 
microscopic theories and models, various generalizations of the Ginzburg-Landau-type 
models to account for high Tc properties such as the anisotropy and the inhomogeneity 
have been proposed and developed. In low Tc superconductors, electrons are thought 
to pair in a form in which the electrons travel together in spherical orbits, but in op- 
posite directions. Such a form of pairing is often called s-wave |llj . However, in high 
Tc superconductors, experiments have produced strong evidence for d-wave pairing 
symmetry in which the electrons travel together in orbits resembling a four-leaf clover 
(cf. [H dni 112 [B] and the references therein) . Recently, the d-wave pairing has gained 
substantial support over s-wave pairing as the mechanism by which high-temperature 
superconductivity might be explained. In generalizing the Ginzburg-Landau models 
to high Tc superconductors, the key idea is to introduce multiple order parameters in 
the Ginzburg-Landau free energy functional. These models, which can also be derived 
from the phenomenological Gorkov equations [B], have built a reasonable basis upon 
which detailed studies of the fine vortex structures in some high T^ materials have 
become possible. We refer the reader to |4] [101 El |6| and the references therein for a 
detailed exposition on modeling and analysis of d-wave superconductors. 

We obtain equation ( |1.1| from the Ginzburg-Landau-type d-wave model consid- 
ered in [4] (also see [lOl I12j ) in absence of applied magnetic field by neglecting the 
zeroth order nonlinear terms but retaining the leading terms. In the equation, u (no- 
tation ij^d is instead used in the cited references) denotes the d-wave order parameter. 
We note that the original order parameter ipd in the Ginzburg-Landau-type model 
[TUl U] is a complex- valued scalar function whose magnitude represents the density of 
superconducting charge carriers, however, to reduce the technicalities and to present 
the ideas, we assume u is a real- valued scalar function in this paper and remark that 
the finite element methods developed in this paper can be easily extended to the com- 
plex case. We also note that the parameter 5 appears in the full model as 5 = — ^, 

where (3 is proportional to the ratio with T^o and Tdo being the critical 

temperatures of the s-wave and d-wave components. Clearly, /3 < (or J > 0) when 
Tso <T < Tdo and /3 \ — oo (or S \0) asT y Tdo- Hence, 6 is expected to be small 
for d-wave like superconductors. 

The primary goal of this paper is to develop conforming finite element methods 
for the reduced d-wave model ( |1.1| . Since the bi-wave term is the leading term in the 
full d-wave model, see [H Section 4], any good numerical method for (1.1) should be 
applicable to the full d-wave model. It is easy to see that the energy space for the 
bi-wave equation (1.1 1 is V {v G H^{fl); □« e L^{fl)} (see Section [2|. Our main 
task then is to construct finite element subspaces V'^ of the energy space V which 
should be as simple as possible but also rich enough to have good approximation 
properties. To this end, we note that H^{il) C V C H^{U), and hence, the desired 
finite element space should satisfy V'^ C V C H^{il.). This immediately implies 
that C C°(ri) (see |31 H]). On the other hand, since ^ is a proper subspace of 
H'^{n), the condition C C°(U) does not guarantee that V'^ C V. Hence, C° 
(Lagrange) finite element spaces are in general not subspaces of V. An intriguing 
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question is what extra conditions are required to make a C*^ finite element space to be 
a subspace of V. To answer this question, on noting that H^{il) C V, one may choose 
such that V'^ C H^{^1), that is, V'^ is a C-^ finite element space such as Argyris 
finite element space (cf. [31 Chapter 6]). Trivially, V'^ C H^{il) C V. It turns out 
(see Section |4]) such a choice would work since it can be shown that the finite element 
solution so defined converges with optimal rate in the energy norm of V. However, 
since finite elements require either the use of fifth or higher order polynomials 
with up to second order derivatives as degrees of freedom [T3J [T3] , or the use of exotic 
elements O Chapter 6] , it is expensive and less efficient to solve the bi-wave equation 



( 1.1 1 using finite elements. This then motivates us to construct low order non-C^ 
finite elements which give genuine subspaces of V and to develop other types of finite 
element methods such as nonconforming and discontinuous Galerkin methods [7]. 

The remainder of the paper is organized as follows. Section [2] contains some 
preliminaries and the functional setting for the bi-wave problem. Well-posedness of 
the problem and regularity estimates of the weak solution are established. Because 
is a hyperbolic operator, the usual regularity shift for fourth order elliptic problems 
does not hold for the bi-wave problem, instead, a weaker shifting "rule" only holds. 
Section [3] devotes to construction and analysis of piecewise polynomial subspaces of 
V. First, we give a characterization of such subspaces. It is proved that a subspace 
of V is "necessarily" a finite element space on a general mesh. However, non-C^ 
finite elements are possible on restricted meshes. Second, we construct two such finite 
elements. The first one is a cubic element and the second is a quartic element. Third, 
we establish the approximation properties for both proposed finite elements. Because 
both elements are not affine families, a technique of using affine relatives (cf. [11 [3]) 
is used to carry out the analysis. Finally, optimal order error estimates in the energy 
norm of V are proved for the finite element approximations of problem ([0J-([L2]) 
using the proposed finite elements. Suboptimal (and optimal in some cases) order 
error estimates in the L^-norm are also derived using a duality argument. In Section 
|4]we present some numerical experiment results to gauge the efficiency of the proposed 
finite element methods and also to validate our theoretical error bounds. 

2. Preliminaries and functional setting. Standard space notation is adopted 
in this paper. We refer the reader to [2l|3] for their exact definitions. In addition, (•, •) 
and (•, ■)gQ are used to denote the L^-inner products on and on d^l, respectively. 
C denotes a generic h and (5-independent positive constant. We also introduce the 
following special space notation: 

Vo := {v e V n H^in); d^v\g^^ = 0}, {v,w)v '-^ 6{Du,Dw) + {v,w), 
\\v\\v ■■= \/{v,v)v- 

It is easy to verify that (•, •)y is an inner product on V, hence, || • ||v is the induced 
norm, and V endowed with this inner product is a Hilbert space. We remark that 



all above claims do not hold in general if the harmonic term Au is dropped in (1.1 1 
because the kernels of the bi-wave operator and the wave operator □ may contain 
non-zero functions satisfying the homogeneous Dirichlet boundary condition [1 



The variational formulation of ([TT|)-([L2| can be derived easily by testing (1.1 1 



against a test function v £ Vq and using integration by parts formulas. Specifically, 
it is defined as seeking u G Vq such that 

(2.1) A'{u,v) = {f,v), 
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where 



A^{u,v) {u,v)v, 

and (•, •) denotes the pairing between V and its dual, V* . 
We now show that problem (2.1 1 is well-posed. 

Theorem 2.1. For any f G V* , there exists a unique solution to (2.1 1. Further- 
more, there holds estimate 



(2.2) 



\u\\v<\\f\\v^. 



Proof. We note for v,w € Vq, 

(2.3) A'{v,v)>\\vry, 

(2.4) \A^{v,w)\ <\\v\\v\\w\\v 

Then, existence and uniqueness follows directly from an application of the Lax- 
Milgram Theorem (cf. [2, 3 ) and using the fact that Vq is a Hilbert space with 
the inner product (•, ■)v- The estimate (2.2) follows from (2.3l and (2.1) after setting 
V — u and w = u. U 

We note iJ^(ri) is a proper subspace of V, so in general u ^ H^{il) if / £ V* . 
However, for smoother function / we have the following regularity results. 

Theorem 2.2. Assume that the boundary of the domain f2 is sufficiently 
smooth. Letsi,S2 be two nonnegative integers. Then there exist constants Cs^^s2t(-'si,s2 
such that the weak solution u of (2.1) satisfies 



(2.5) 

V~5\\U^dl}dl 
(2.6) 



d^d^fW- 



^~5\\VUdl}dl-u\\L^ 



ifdl}dl^f&V* 



ifdi}di^f^L\n). 



> 



Proof. First, we consider the case that u and / have compact support. Let w := 
d^^dy^u and g := d^^dy^f. Because equation ( |1.1| is a linear equation, differentiating 
the equation immediately verifies that w and g satisfy 



(2.7) 



Sn'^w - Aw 



that is, w is a solution of the bi-wave equation with the source term g. Since u is 
assumed to have a compact support, then w also satisfies the homogeneous boundary 
conditions in (1.2). Thus, it follows from Theorem 2.1 that 

\\w\\v < hWv, 



which gives (2.5 1 with Cg^^s^ — 1. 

To show (2.6), it suffices to prove that 

(2.8) VS\\Vnw\\L2 + \\Aw\\l2 < Cs„sM\l- 

S2 



which is equivalent to prove that (2.6) holds for si = S2 = 0. To this end, testing 
-S(D^w, Aw) + \\Aw\\l2 = -{g, Aw). 



(2.7) with — Aw yields 



Using the following integral identity 



followed by using Green's identity (for A) in the first term on the left hand side we 
get 

-(5(n2w, Aw) = 5\\VUw\\l^. 

Here, we have dropped the boundary integral terms because w has a compact support. 
Combining the above two identities for w and using Schwarz inequality yield 



5\\vnw\ 



L2 



l|Au;||i. <^||ff||i. 



JllA^lli.. 



Hence, the above inequahty and (2.7 1 imply that (2.8 1 holds with Csi.s2 — 2\/2+ 1 



Second, in the case u and / do not have compact support, it is clear that w and 



g still satisfy (2.7 1. However, w and its derivatives may not satisfy the homogeneous 



boundary conditions in (1.2 1. To get around this difficulty, the well-known tricks are 



to use the cutoff function technique (see [SJ [5] ) for interior estimates and to use the 
flattening boundary technique for boundary estimates. The cutoff function technique 
involves testing (2.7) by wS, and — Aw^, instead of w and — Aw, for a smooth cutoff 



function ^. Integrating by parts on the left hand side and using Schwarz inequality 
and the properties of the cutoff function then yield the desired interior estimate similar 
to (2.5 1 and (2.6). The flattening boundary technique involves locally mapping the 



curved boundary into a flat boundary by a smooth map (this requires the smoothness 
of the boundary d^). After the desired boundary estimates are obtained in the new 
coordinates, they are then transferred to the solution w in the original coordinates. 
We omit the technical derivations and refer the interested reader to |S][5] for a detailed 
exposition of these techniques applying to other linear PDEs. □ 

3. Construction and analysis of finite element methods. 

3.1. Characterization of finite element subspaces of V . Let 7/i be a quasi- 
uniform triangulation of Q. with mesh size h G (0,1), and for a fixed T e 7^, let 
(A^, AJ, A^) denote the barycentric coordinates, and (1 < i < 3) denote the vertices 
of T. We also let (1 < ? < 3) denote the edge of T of which a,j is not a vertex, and 
bi denote the midpoint of edge e^. Define the interior and boundary edge sets of Th 



We also set 



£kyj£^ 



h I 



and for T (^Th, 



uj{T) closure ( |J T' 

For any e e such that e = Ti n T2, and v e H^{Ti)r] H^{T2), define the jumps of 
V across e as (assuming the global label of Ti is bigger than that of T2) 



„Ti| 



where v^^ = v\rp , and :— v^^]^ if e e . 

Similarly, for v G H^(Ti)n H^{T2), a £ R?, we define the jumps of daV :— Vv ■ a 
as follows: 

[do.v]\^ := dc.v'^'l^ - da^v^'l e = dnn dT^ € £{, 
We also define the shorthand notation 



In the rest of the paper, we shall often encounter the following characterization 
of the meshes. 

Definition 3.1. For e £ £h, let n and t denote the outward unit normal and 
unit tangent vector of e, respectively. We say that e is a type I edge if 

(3.1) n = T or n — ~T. 



Otherwise, e is called a type II edge if condition (3.11 does not hold. 

Remark 3.1. (a) If e is a type I edge, then n = (rii,— rt2) = ±t = ±(Ti,r2) = 
^{'^'2, —ni). Therefore, 

(±1,±1). 

That is, the edge e makes an angle of ^ in the plane with respect to the x-axis. 
Examples of meshes such that every triangle in the partition has exactly zero and one 
type I edges are shown in Figure \371\ and examples of meshes such that every triangle 
has exactly two type I edges are shown in Figure \3.S\ 

(h) For T G 7^, gj C dT, let rS^^ and r^*' denote the outward (from T) unit normal 
and unit tangent vector of Ci, respectively. Then using the formula 



we conclude that Ci is a type I edge if and only if 

VXT ■ VA^ = 0. 



To construct finite element subspaces of V^, we first provide the following two 
lemmas, which characterize such spaces. 

Lemma 3.2. Let he a suhspace ofV consisting of piecewise polynomials, and 
suppose there exists a type II edge e G with e = dTi n dT2. Then for v e X'^ , there 
holds the inclusion v G H^(Ti U T2). 

Proof. Since X^ is finite-dimensional with X'^ C II^{il), we have the inclusion 
X'^ C C°(n). We also note that it suSices to show v e C^(Ti U T2) for any v e X^, 
which in turn is equivalent to show 

Let n and r denote the normal and tangential direction of e, respectively. Rewrit- 
ing V as 



V^{ve H\n); Vu e i7(div; ft)}, 
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Fig. 3.1. Example of meshes of the domain Q — (0, 1)'^ such that every triangle has no type I 
edges (left), and one type I edge (right). 




there holds for v ^ X^^ 

Next, using the assumption n ^ ±t, we can write for any constant vector a G 



[dai 



1 — (r • n)' 



But [drv]\^ = since v G C°(ri) and is a polynomial of one variable. Hence, 



[9. 



implies that [dav] 



= 0. The proof is complete. □ 
Corollary 3.3. Suppose is a subspace of V consisting of piecewise polyno- 
mials, and suppose there exists no type I edges in the set Then C H^{^). 

Lemma 3.4. Suppose Et is a linearly independent set of parameters uniquely 
determining a kth-degree polynomial v on an interior triangle T € that includes 
only function and derivative degrees of freedom. Suppose further that v is continuous 
in io{T), □« G L'^{uj{T)), and T has at least two type II edges that are in the set E^. 
Then fc > 5. 

Proof. If T has three type II edges, then by Lemma 3.2 v G II'^{uj{T)), and it 
follows that fc > 5 (cf. [3, p. 108], also see p^fli]). 
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Suppose T has exactly two type II edges, without loss of generality, assume ei is 
type I. By the proof of Lemma 3.2 v is across edges 62 and 63. Let denote the 
order of prescribed derivatives at vertex in the set Ey, let denote the number 
of function value (or equivalent) degrees of freedom in the set on edge Ci, and let 
Si denote the number of (non-tangential) directional derivative value (or equivalent) 
degrees of freedom in the set St on edge e^. Since v is continuous in lo(T), we have 

(3.2) M2 + M3 + ™i > - 1, 

Ml + M3 + ^2 > fc - 1, 
Ml + A*2 + ^3 > - 1, 

and since is continuous across 62 and 63, 

(3.3) Mi+M2 + S3>fc, 

Ml + M3 + S2 > fc. 

Adding up the above five inequalities yields 

4/ii + 3/i2 + 3/i3 + mi + 7712 + TO3 + Si + S2 > 5fc — 3. 

Because the set is linearly independent, and the dimension of equals 
(^±i^, there holds 

(3.4) ik+l){k + 2) ^ ^ 1 1 ^ ^^^^^ ^^^^ ^^^^^^ 

i=l 
3 1 

- 51 2 ''^^ ^ + 2) + S/c - 3 " 4/ii - 3/i2 - 3/i3- 

i=l 

Thus, 

(3.5) (fc2 - 7fc + 8) > (Ai? - 5Aii + 2) + (^^2 - 2)(m2 - 1) + (M3 - 2)(m3 - !)• 

It is clear that k must be greater than two, therefore, it suffices to show that k 
cannot equal three or four. 



Case fc = 3; If fc = 3, by ( |3.5| we get 

(//I - - 2) + {^l2 - 2)(m2 - 1) + (m3 - 2)(m3 - 1) < 0, 

and since /i^ are integer- valued, we have 

1 < M3 < 2, 1 < M2 < 2, 2 < ^1 < 3. 



But by (3.4 1, we immediately obtain 

3 



(fc + l)(A:-t-2) ^l, 



10 . 2 

i=i 



which is a contradiction. 

Case k — 4; As in the previous case, if fc = 4 we have 

ifii - 3)(Aii - 2) + (Ai2 - 2)(m2 - 1) + (M3 - 2)(m3 - 1) < 0. 



Since 



1 < Ai3 < 2, 1 < Ai2 < 2, 2 < /ii < 3, 

and 



15 



(/c + l)(fc + 2) 



2 - ^ 2 

1=1 



3 1 



it is not hard to check that there can only be the following three subcases: 
(3.6) (^i,/i2,/i3) = (2, 1,2), (^i,^2,A*3) = (2,2, 1), (^i, ^2, A^s) = (2, 1, 1). 

If the first subcase holds, then all degrees of freedom lie on the vertices, therefore. 



0, 1 < i < 3. However, it follows from (3.3 1 that 

3 = /xi + ^2 > 4, 

which is a contradiction. 

A similar argument can be used to exclude the second subcase in (3.6 1. Now, 
suppose /zi = 2, — 1, and — 1. By (3.2) and (3.3l, we have 

W3 > 1, S3 > 2, S2 > 1. 

But this implies that 

3 1 

15 > ^ {-(/i, + + 2) + im] + S2 + S3 > 16, 

i=l 

a contradiction. Thus, the third subcase can not happen, either. Therefore, we must 
have k > 5. The proof is complete. □ 

By Lemmas |3.2| and |3.4[ and Corollary |3.3| we conclude that unless certain types 
of meshes are used, we must resort to either finite elements such as Argyris, Hsieh- 
Clough- Tocher, Bogner-Fox-Schmit elements (cf. [HE]), or special exotic elements 
(e.g. macro elements), or nonconforming elements (cf. |3) to solve problem (1.1)- 
( |1.2[ ), However for special meshes, we now show in the following subsections that it is 
feasible to construct low order finite element subspaces of V. 

3.2. A cubic conforming finite element. To construct a cubic conforming 
finite element, we assume that T/j is a triangulation of il and every triangle of 7/j has 
two type I edges. Examples of such meshes are shown on a square domain in Figure 



3.2 Our cubic finite element 6*3 := {T, Pt,T,t) is defined as follows: 

(i) T is a triangle with two type I edges, 

(ii) Pt = ^3(1"), the space of cubic polynomials on T, 
v{ai) 1 < * < 3, 
viais) 1 < i < 2, 
Vti(ai) • (flj — Gi) I < i <2, 1 < j < 3, j 7^ i, 

where 63 is a type II edge. 
Lemma 3.5. The set St is unisolvent. That is, any polynomial of degree three is 
uniquely determined by the degrees of freedom in Y.t- 
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(iii) Si 



Fig. 3.3. Element . Solid dots indieate function evaluation, circles indicate first derivative 
evaluation, and arrows indicate evaluation of derivatives in the direction n. 



Proof. Suppose V E PsiT) equals zero at all the degrees of freedom in St- To 
complete the proof, it suffices to show v = since dim(P3(T)) ~ dim(ST) = 10. 

Recall that 63 is a type II edge, ei and 62 are type I edges of T. Let Wi be the 
restriction of w on C dT as a function of a single variable, then Wi is a polynomial 
of degree three which satisfies 

u;:(0) = u;,(0) = =u;,(l) = i=l,2, 

In either case, we conclude Wi = 0. 

Next, let Z3 be the restriction of dnV on 63 as a function of a single variable. Then 
Z3 is a polynomial of degree two satisfying 

^3(0) = ^3(^) = ^3(1), 

which then infers Z3 = 0. 

From the above calculations, we conclude that (A^")^, , and A^ are factors of 
V. However, this is not possible, since u is a polynomial of degree three, unless v = 0. 
The proof is complete. □ 

Let be the finite element space associated with 6*3 , that is, 

^3 — G P3(r), V is continuous at every degree of freedom in St, VT e T^}. 

We now show that V^^ is a subspace of V. 

Theorem 3.6. There holds the inclusion Vt_ C V. 



Proof. Let v e V^. By the proof of Lemma 3.2 it suffices to show v and dnV are 
both continuous across interior edges of 7^. Let Ti and T2 be two adjacent triangles 
with common edge e, and w be the restriction of v"'"^ — v^^ along e as a function of a 
single variable. We then have 

w'(0) wi(0) = w(^) = ^(1) = if e is type I, 
w'{Q) = w(0) = w(l) = w'{\) = if e is type II. 



Thus, w = and the inclusion V^' C C°(r2) C H^{n) holds. 
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Next, we observe that if e is a type I edge, then 



[dnv]\^=±[drv]\^^=0. 



Hence, dnV is continuous across e. On the other hand, if e is a type II edge, let z be 
the restriction of [9fjt;]|e = dnv^^ — dnv'^^ along e as a function of a single variable. 
Since 



z(0) 



4> 



and z is a polynomial of degree two. it follow that [3fiw]|^ = 0. So dnV is also 
continuous across e. This then concludes the proof. □ 

Remark 3.7. We note that V^" <t H'^{0.) because t C^(n). 

3.3. A quartic conforming finite element. In this subsection, we again as- 
sume that 7^ is a triangulation of and every triangle of Th has two type I edges. 
We then define the following quartic finite element 5*4 := (T, Qt,'^t)'- 

(i) T is a triangle with two type I edges, 

(ii) Qt = Vi{T), the space of quartic polynomials on T, 
v{ai) 1 < « < 3, 
"(ftjjg), ■(;(ai33) 1<«<2, 

w(ai23), 

9ftW(aii2), dnV{ai22), 

where ea is a type II edge, and aijt, = | (oj 



(iii) S3 



1 < i < 3, i ^ i, 
J + at). 




Fig. 3.4. Element . Solid dots indicate function evaluation, circles indicate first derivative 
evaluation, and arrows indicate evaluation of derivatives in the direction n. 

Lemma 3.8. The set St is unisolvent. That is, any polynomial of degree four is 
uniquely determined by the degrees of freedom in Sy. 

Proof. Suppose v E P4(T) equals zero at all the degrees of freedom in S^, and let 
Wi be the restriction of z; to Cj as a function of a single variable. Then 

w'iiO) = Wi{0) = Wi{^) = Wi{^) = Wi{l) = i = 1, 2, 

w'M = W3{0) = wsi^-) = «;3(1) = w',{l) = 0. 
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Thus, Wi = 0,i — 1,2,3. 

Next, letting be the restriction of dnV on 63 as a function of a single variable, 
we have 

Z3{0) = z,{^) = z,{'^) = z,{l) = 0. 

Hence, Z3 = 0. 

From the above calculations, we conclude that v = aAiA2A3 for some a G R. 
However, since = ■i;(ai23) = f^, we have a = 0. The proof is complete. □ 

Theorem 3.9. Let he the finite element space associated with S'4 , that is, 

V4 ~ {v\t e P4(T), V is continuous at every degree of freedom in St, VT G 7/i}. 

Then there holds the inclusion C V . 

Proof. Let v G , and suppose T^,T'^ G 7^ are two adjacent triangles with 
common edge e. Let w be the restriction of [v]\e = v^'^ — v^^ along e as a function of 
a single variable, from 

1 2 

w'(0) = u;(0) = w{-) = w{-) = w(l) ==0 if e is type I, 

o o 

w'{{)) = w(0) = u)(^) = w{l) = w'(l) = if e is type H, 

we conclude w = Q. Hence, the inclusion C C°(ri) C iJ^(ri) holds. 

If e is a type H edge, we let z denote the restriction of [9fit;]|e = dnv^'^ — dnv'^^ 
along e as a function of one variable. It follows from 

^(0) = ^(^) = ^(^) = ^(1) = 0, 

that 2; = 0. 

Finally, if e is a type I edge, we use the fact that v is continuous to conclude 

[dnv]\^ = ±[drv]\^ = 0. 

Thus, Vi" CV.D 

Remark 3.10. We note that ^ H^{fl) because Vl' (t C^fl). 

3.4. Approximation properties of the proposed finite elements. Let 

H^u G Pfc(T) denote the standard interpolation of v associated with the finite ele- 
ment S^, and define H^w G such that H^vj^, = H^(t)|^), VT G Before stating 
the approximation properties of the interpolation operator H^, we first establish the 
following technical lemma concerning the mesh 7^. 

Lemma 3.11. Suppose T G 7^ has two type I edges, and without loss of generality, 
assume 63 C dT is a type II edge. Then there exists a constant C > that depends 
only on the minimum angle ofT such that 

where (3^^^ = r^^^ • n^^^ 

Proof. Since both type I edges of T make an angle of ^ with respect to the a;-axis 
(cf. Remark 3.1 ), then there exists G (0, ^] such that the angles of T are ^, 9, and 

l-e. 
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X 



Fig. 3.5. Embedding T into an isosceles triangle 



Next, we embed T into an isosceles triangle as shown in Figure [375] and then 
obtain 

^(3)^(M, r(3).n(3)^z2^, x^cosC--6)z, y = sin(--0)z. 
z 4 4 

Hence, 

^(3) . ^(3) ^ ^ -2sin(- - 9) cos(- - 0) = -cos(20), 

z^ 4 4 

which implies that 

1 - (/3(3))2 ^ 1 _ . ri(3))2 = 1 _ cos2(2^^) = sin2(20). 

The proof is complete. □ 

Remark 3.12. IfTh is a uniform criss-cross triangulation offl, then l — {(3^^^)'^ = 
1 for all type II edges, 63. 

The next theorem establishes the approximation properties of the proposed cubic 
and quartic finite elements. 

Theorem 3.13. For all m > 0, p, q € [l,oo] which are compatible with the 
inclusion 

there holds 

(3.7) ^-ni^«||w"".'.(T) <c4^'""^'"^||^;||v^.+i.p(T) Vw G T^'^+i'f (T), 

where Ht = diam(T). 

Proof. The case 5*3; Since is not an affine family in general, the standard 
scaling technique can not be used directly to prove (3.7). To get around this difficulty, 
the trick is to introduce an affine "relative" of 6*3 and to estimate the discrepancy 
between S'3 and its "relative" . To this end, we introduce the following element 53 := 

(r,PT,S^): 

(i) T is a triangle with two type I edges, 

(ii) Pt = P3(r), 

v{ai) 1 < i < 3, 

v{aik) I <i <2, 

Vv{ai) ■ (flj - ftj) 1 < z < 2, 1 < j < 3, j 7^ *7 
Wv{b3) ■ (03 - 63), 
where edge 63 is of type II. 
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(iii) 




It is easy to see that is unisolvent in P3(r), and that any two triangles are 
affine equivalent. Therefore for allp, g e [1, oo], < m < 4 with W^'P^T) ^ W^™'«(r), 
there holds [3] 

(3.8) \\v-A^v\\w^.,(^T)<Chl'"'^~''^\\v\\w^.p^T) yveW^-P{T), 

where is the interpolation operator associated with ^3. 

Define := U'^ - , and note that for v e W^''p{T), &Jv\^ = for i = 1,2,3. 
Consequently, 

Vvih) ■ (as - 63) = ^—^^[{a, - 63) • {n^'^ - r^^)/?^^))^,,^, {v - Ajv){bs)}, 

where n^^) = {nf\-n'i^), (i^^) t(3) . nO), ^(3) = and r^^) denote 

respectively the unit normal and tangential direction of edge 63. 

Next, let (73 be the basis function associated with the degree of freedom Vu(63)(a3 — 
63) in Ti'j,. We then have 

er.= ^3^{(a3-63) . (n(3) -r(3);3(3))a,<3,(.~Aj.)(63)}g3. 
Therefore, 

l|e^Hlw™..(T) < - • 1"^'^ - • \\v ~ AlvU..^^r)\\<lz\\w^:ij) 



Finally, by ( |3.8| and Lemma 3.11 we get 

l-(/3(3))2>c, 103-63! <C/ir, 

1^(3) _ ^(3)^(3) I < 2, 11^;- A^«||v^.,o.(^) < c4"^||^;|1h.4,p(t), 

lk3|lw".'J(T) < C'/iT ' ' 

where C only depends on the minimum angle of T . Hence, 



2 2 
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and consequently, 

\\v - Il^v\\w^-.i{T) < \\v - A3 w|liv™.<!(r) + ll0l'w|iiV"'>9(T) < Ch^j, " \\v\\wi.p{T)- 



The case S^: We use a similar argument to show (3.7) for the element S!l. First, 
we introduce the following "relative" S'4 := {T,Qt,'^t) of 



(i) 
(ii) 



(iii) 



T is a triangle with two type I edges, 

Qt = P4(r), 

v{ai) 
v{ai23), 

Vv{ai) ■ {a.j - ai) 
Vw(aii2) • (03 - 0112), 
V-y(ai22) • (03 - 0122), 
where edge 63 is of type II. 



1 < i < 3, 
1 <i<2, 



1 < « < 2, 1 < j < 3, 




edge e^ 



Fig. 3.7. Element S'^. 

Next, let A4 be the interpolation operator associated with S'^, and set 64 := 
nj — A4". Let ri be the basis function of the element ^4 that is associated with 
the degree of freedom Vi'(aii2)(a3 — 0112), and let r2 be the basis function that is 
associated with the degree of freedom Vt;(ai22)(a3 — 0122). Then for v E W^'P{T) 

« = l_()j(3))2 {(«3 - aii2) • {n^'^ - T(3)/3(3))a^(3, {v - Alv)ian2)ri 

+ {as - 0122) • {n^^^ - T(3)/?(3))a^,3, {v - Aft-) (0122)^2}. 



Using the fact ^4 is affine equivalent and applying Lemma 3.11 we get 

1 - > c, |a3 - aii2|, |a3 - 0122! < Chr, 

1^(3) _ ^(3)^(3) I < 2^ \\v~Alv\\w^^^^T) < c4~^lkllw=..(T), 

— m-{- - 

WriWw^.iiT) < Chj, ',1=1,2. 



Therefore, 



\Qi\\w'^'iiT) < Chrp " " ||w||vi/!3.p(T)> 
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and consequently, 

\\v - n![||H"".9(T) < \\v - Alv\\w^.i{T) + ll©! "II W".?(T) < Ch^. '"^^ " ||w||w5.P(T)- 

The proof is complete. □ 

We note that if a uniform criss-cross mesh is used such that every triangle has 



two type I edges (see Figure 3.2), then Vw(63)(a3 — 63) = ±dnV in the definition of 
T,'rp. This observation leads to the following corollary. 

Corollary 3.14. Suppose is the uniform criss-cross triangulation offl, then 
53 = 53 . Hence, S!} is an affine family. 

4. Finite element formulation and convergence analysis. Let V^'* (fc = 
3, 4) be the finite element subspaces of V constructed in the previous section. Define 

Based on the weak formulation ( |2.1| , we define our finite element method for problem 
(1.1)-(1.2| as seeking Uh & V^g such that 



(4.1) A'{uh,Vh)^{f,Vh) yvheV^o- 



On noting (2.3l-(2.4), an application of Cea's Lemma |3J yields the following 
result. 



Lemma 4.1. There exists a unique solution to (4.11. Furthermore, the following 
error estimate holds: 

\\u-Uh\\v<C inf \\u-v,i\\v. 

Combining Lemma |4.1| and Theorem |3 . 1 3| with p = q = 2, m = 1,2 we immedi- 
ately get the following energy norm error estimate. 
Theorem 4.1. Ifu e H^in) (s > 3) then 

\\u-Uh\\v < Ch'^-^{Vd + h)\\u\\Hi, ^ = min{fc + l,s}. 



Next, using a duality argument, we obtain an error estimate in the L^-norm. 
Theorem 4.2. Suppose u G H'^(fl) (s > 3). Then there holds the following error 
estimate: 

(4.2) Wu-Uhh^ <CCoflh'^-\Vd + h)\\u\\H'- ^ = min{A: + 1, s}. 



Proof. Denote the error by := u — u^, and let ip Vq he the solution to the 
following auxiliary problem: 

A^iip,v) = {eh,v) yveVo. 



It follows from Theorems |2 . 1 1 and |2. 2| that the above problem has a unique solution ip 
and 

(4.3) VS\\Vn^\\L2 + \\A^\\l- < C'o.ollehllL^. 
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We then have 

(4.4) \\e,,\\l,^A\eh,v) = A\eH,^-Vl:v) < WeM^-VH^Wv, 

where Vj^ denotes the L^-projection to V^q. 
By the definition of || • \\v and (|4.3| we get 



(4.5) ||<^ - viicpWv < V~s\\Dip - vtn^u^ + wv^ - p.^v^iu. 

< C\/5h\\\IULp\\L2 + C/i||V(V(^)||l2 

< C/i(V(5||Vn^||i2 + ||A^||i2) 

< CC(3flh\\eh\\L^. 



Thus, it follows from Theorem |4.1[ ( |4.4[ ), and ( |4.5| that 

||e,.|U2 < CC•o,o/l'"'(^+/l)lklk- 
The proof is complete. □ 

We conclude this section with a few remarks. 

Remark 4.3. (a) The energy norm error estimate is optimal, on the other hand, 
the and norm estimates are optimal provided that VS ~ h. 

(h) All above convergence results only hold for the restricted meshes, that is, every 
triangle of the mesh Th needs to have two type I edges. As already mentioned at the 



end of Section 3.1 for arbitrary mesh Th, C V will implies that (and V^) 
needs to be a finite element space on Th such as Argyris, Hsieh-Clough-Tocher, 
Bogner-Fox-Schmit elements (cf. JEI)- In such a case, it follows from Lemma \4-. l\ that 

\\u-Uh\\v<C inf \\u-Vh\\v 

< C inf {^/5\\u-Vh\\m + \\u-Vh\\m] 

< Ch'^'\V6+h)\\u\\He, 

where £ — min{fc + l,s} and k{> 5) is the order of the finite element. Thus, we 
still get optimal order error estimate in the energy norm. Although, as expected, using 
finite elements is not efficient to solve the bi-wave problem (cf. 171). 

5. Numerical experiments and rates of convergence. In this section, we 
provide some numerical experiments to gauge the efficiency and validate the theoret- 
ical error bounds for the finite element S'3 developed in the previous sections. 

Test 1. For this test, we calculate the rate of convergence of ||m — for fixed 5 
in various norms and compare each computed rate with its theoretical estimate. All 
our computations are done on the square domain — (0, 1)'^ using the criss-cross 
mesh. We use the source function 

f{x,y) =- 20487r'*5(cos^(47ra;) - sm^{4:7:y)) - 327r2| sin^(47ry)( cos^(47rx) — sin^(47ra;)) 
+ sin^(47ra;)(cos^(47ry) - sin^(47r?;)) |, 

so that the exact solution is given by u{x, y) — sin^(47ra;) sin^(47ry). 
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We list the computed errors in Table O for (5-values 10,1,10^2 and 10"^ and 
also plot the results in Figure [5T2| As expecte d, th e rates of convergence depend on 



4.1 



tells us that for v (5 > > /i 



both the parameter h and 6. In fact, Corollary 

\\u - UhWv < Ch^{Vd + h)\\u\\Hi < Ch^WuWHi, 
\\u - UhWm < Ch^VS + h)\\u\\H^ < Ch^\\u\\m, 
\\u ~ u/i||l2 < CCo.oh^iVS + /i)||m||h4 < CCo,oh^\\u\\H*, 



while for V S < h 

\\u - UhWv < Ch^iVS + h)\\u\\H^ < Ch^Mn^, 
\\u - Whilifi < Ch^iVS + h)\\u\\Hi < Ch^\\u\\Hi, 
\\u~Uh\\L2 <CCo,oh''{V6 + h)\\u\\H^ <CCo,oh^u\\H^. 

We find that the computed bounds agree with these theoretical bounds. 

In addition, although a theoretical proof of the following convergence rate has yet 
to be shown, the computed solutions also indicate that 



where 



ll"" - Uh\\2ji < Ch{\/5 + h)\\u\\Hi, 




Fig. 5.1. Test 1. Computed solution (left) and error (right) with (5—10 ^ and h — 0.01. 

Test 2. This test is the same as the first, but we now use the following source 
function: 

/ = !■ 

We note that the exact solution is unkn own. We plot the solution with h = 0.01 and 
<5-values 10,1,10~^, and 10"^ in Figure 5.3 As expected, the solution is more and 
more like the solution of the corresponding Poisson problem as 5 gets smaller and 
smaller. 
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S h 


• 11^2 err. (cnv. rate) 


• Wjji err. (cnv. rate) 


• II ft, 2 err. (cnv. rate) 


• \\v err. (cnv. rate) 


10 0.5000 


4.17(-) 


26. 4(-) 


311.62(-) 


2191. 62(-) 


0.3333 


2.76E-01(6.694) 


9.54(2.514) 


284.05(0.228) 


1147.49(1.596) 


0.2000 


1.59E-01(1.079) 


2.99(2.273) 


211.50(0.577) 


535.69(1.491) 


0.1000 


4.41E-03(5.176) 


3.75E-01(2.995) 


52.54(2.009) 


153.70(1.801) 


0.0500 


4.64E-04(3.248) 


9.11E-02(2.041) 


21.61(1.282) 


39.84(1.948) 


0.0400 


2.31E-04(3.117) 


5.82E-02(2.010) 


16.79 (1.130) 


25.61(1.980) 


0.0200 


2.79E-05(3.054) 


1.45E-02(2.004) 


8.05(1.060) 


6.44(1.992) 


0.0100 


3.45E-06(3.014) 


3.62E-03(2.001) 


3.98(1.016) 


1.61(1.998) 


0.0083 


1.99E-06(3.006) 


2.52E-03(2.000) 


3.31(1.006) 


1.12(1.999) 


0.0067 


1.02E-06(3.004) 


1.61E-03(2.000) 


2.65(1.004) 


0.72(1.999) 


1 0.5000 


3.93(-) 


25. 3(-) 


306. 88(-) 


238. 43(-) 


0.2500 


2.75E-01(3.837) 


9.52(1.413) 


283.58(0.114) 


123.22(0.952) 


0.2000 


1.57E-01(2.523) 


2.98(5.210) 


210.95(1.326) 


56.24(3.515) 


0.1000 


4.40E-03(5.152) 


3.75E-01(2.989) 


52.53(2.006) 


15.71(1.840) 


0.0500 


4.63E-04(3.249) 


9.09E-02(2.043) 


21.58(1.284) 


4.07(1.950) 


0.0400 


2.31E-04(3.118) 


5.81E-02(2.012) 


16.76(1.132) 


2.61(1.981) 


0.0200 


2.78E-05(3.055) 


1.45E-02(2.005) 


8.03(1.061) 


0.66(1.992) 


0.0100 


3.44E-06(3.015) 


3.61E-03(2.001) 


3.97(1.016) 


0.16(1.998) 


0.0083 


1.99E-06(3.006) 


2.51E-03(2.000) 


3.31(1.006) 


0.11(1.999) 


0.0067 


1.02E-06(3.004) 


1.61E-03(2.000) 


2.64(1.004) 


0.07(1.999) 


0.0056 


5.89E-07(2.999) 


1.12E-03(2.000) 


2.20(1.003) 


0.05(1.998) 


10~^ 0.5000 


2.15(-) 


15.4(-) 


276. 56(-) 


17.4(-) 


0.3333 


2.25E-01(3.259) 


8.38(0.879) 


260.32(0.087) 


9.48(0.877) 


0.2000 


1.02E-01(3.556) 


2.53(5.365) 


183.36(1.571) 


3.07(5.047) 


0.1000 


4.21E-03(4.597) 


3.69E-01(2.780) 


52.08(1.816) 


5.22E-01(2.558) 


0.0500 


4.36E-04(3.269) 


8.55E-02(2.107) 


20.31(1.358) 


1.25E-01(2.058) 


0.0400 


2.15E-04(3.175) 


5.38E-02(2.082) 


15.52(1.207) 


7.93E-02(2.049) 


0.0200 


2.50E-05(3.101) 


1.30E-02(2.049) 


7.21(1.106) 


1.94E-02(2.030) 


0.0100 


3.06E-06(3.033) 


3.21E-03(2.016) 


3.53(1.030) 


4.82E-03(2.010) 


0.0083 


1.77E-06(3.013) 


2.23E-03(2.006) 


2.94(1.012) 


3.35E-03(2.004) 


0.0067 


9.02E-07(3.009) 


1.42E-03(2.005) 


2.34(1.008) 


2.14E-03(2.003) 


10-^ 0.5000 


3.93(-) 


23.3(-) 


374.18(-) 


23. 3(-) 


0.2500 


2.28E-01(4.108) 


7.25(1.686) 


233.21(0.682) 


7.25(1.686) 


0.2000 


9.68E-02(3.831) 


1.93(5.929) 


149.75(1.985) 


1.93(5.929) 


0.1000 


3.70E-03(4.708) 


2.81E-01(2.782) 


45.13(1.731) 


2.81E-01(2.782) 


0.0500 


4.92E-04(2.914) 


5.21E-02(2.429) 


13.09(1.786) 


5.21E-02(2.429) 


0.0400 


2.35E-04(3.298) 


2.89E-02(2.647) 


8.54(1.915) 


2.89E-02(2.647) 


0.0200 


1.91E-05(3.626) 


4.11E-03(2.813) 


2.15(1.989) 


4.11E-03(2.813) 


0.0100 


1.28E-06(3.898) 


5.39E-04(2.931) 


0.53(2.024) 


5.39E-04(2.931) 


0.0083 


6.19E-07(3.981) 


3.15E-04(2.943) 


0.36(2.035) 


3.15E-04(2.943) 


0.0067 


2.53E-07(4.005) 


1.64E-04(2.934) 


0.23(2.039) 


1.64E-04(2.933) 



Table 5.1 

Test 1. Errors with estimated rates of convergence 
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